A gauge invariant cluster algorithm for the Ising spin glass 
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The frustrated Ising model in two dimensions is revisited. The frustration is quantified in terms 
of the number of non-trivial plaquettes which is invariant under the Nishimori gauge symmetry. 
The exact ground state energy is calculated using Edmond's algorithm. A novel cluster algorithm 
is designed which treats gauge equivalent spin glasses on equal footing and allows for efficient 
simulations near criticality. As a first application, the specific heat near criticality is investigated. 
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Spin glasses |1| are magnetic materials in which the 
magnetic moments are subject to ferromagnetic or anti- 
ferromagnetic interactions, depending on the position of 
the moments within the sample. The system is frustrated 
in the sense that the arrangement of spins which mini- 
mizes the total energy cannot be determined by consid- 
ering a local set of spins. Stated differently, the change 
of a single spin might cause a reordering of many spins 
when the system relaxes towards a new minimum of en- 
ergy 0- Spin glasses undergo a freezing transition to a 
state where the order is represented by clusters of spins 
with mixed orientations. The relaxation times towards 
equilibrium are typically very large, which impedes effi- 
cient simulations. 

Many efforts have been undertaken to explore equilib- 
rium properties of spin glasses by means of Monte Carlo 
simulations JLiSsil J^J- Thereby, many insights 

have been obtained from the most simple case of the 2d 
Ising model on a square lattice. For the discrete model, 
the bond interactions take values Ji = ±1 at random, 
and the model is characterized by the the probability p 
of finding an anti- ferromagnetic interaction, Jg = —1, at 
a given bond. At present, the existence of a spin glass 
transition in these models at 7^ is still under debate. 

As first noticed b y B ieche et al. and further elabo- 
rated by Nishimori [l^, [131 1 the Ising model with random 
distribution of anti-ferromagnetic bonds has a hidden Z2 
gauge symmetry. As discussed below, this symmetry im- 
plies that gauge invariant observables such as the thermal 
energy or the specific heat are unchanged by a certain re- 
distribution of the anti-ferromagnetic bonds (which may 
also change their number considerably). By exploiting 
this invariance, Nishimori was able to obtain exact results 
for special values of the parameters p and T j3 ^3 ■ 

The generic difficulty in simulating spin systems is 
that the auto-correlation time r increases rapidly with 
the physical correlation length ^ of the system, r cx , 
where z is the dynamical critical exponent. For all local 
update algorithms, z is as large as 2. This is particu- 



larly problematic at small temperatures, when ^ reaches 
the extension L of the lattice and the generation of in- 
dependent configurations in a Markov chain becomes 
extremely cumbersome. As a consquence, the auto- 
correlation times must be monitored very carefully or the 
algorithm might fail to be ergodic. For a pure ferromag- 
net, the ground state is known explicitly: A state with all 
spins parallel minimizes the energy. This knowledge of 
the true ground state can be used to design an efficient al- 
gorithm which microcanonically changes clusters of spins 
with the same orientation. Indeed, the so-called cluster 
algorithms [T^ ITU largely alleviate the auto-correlation 
problem: the dynamical critical exponent drops to val- 
ues z w 0.4 which renders practical simulations on large 
lattices feasible. 

So far, cluster algorithms for the frustrated Ising model 
do not exist. This is mainly due to the fact that the true 
ground state (and hence the structure of the physical 
clusters) is unknown for a generic distribution of anti- 
ferromagnetic bonds. In fact, finding the ground state 
is an NP hard problem in d > 3. For the special case 
d — 2, Edmond's algorithm [IgL IT^ provides a method 
which computes the exact ground state in polynomial 
time. This can be used to clearify the structure of the 
physical clusters in special gauges, which would otherwise 
be obfuscated by the hidden gauge symmetry. 

In this letter, we quantify the amount of frustration in 
the 2d Ising model in a gauge invariant way by count- 
ing the fraction p of vortices (non-trivial plaquettes) in 
a given bond distribution. First, we determine the exact 
ground state energy as a function of p using Edmond's 
algorithm. In order to treat the system near criticality, 
we present a novel cluster update algorithm which pro- 
poses clusters in a gauge independent way. The specific 
heat as a function of the inverse temperature is explic- 
itly evaluated for models with different frustrations, and 
gauge independence is verified. 

The partition function of the frustrated Ising model 
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involves a summation over all spin configurations {cx} 



{<Tx} 
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where the spins located at the sites of the lattice take 
the values ax — ±1- The sum in the exponent extends 
over all bonds t — (xy) and the coupling constants (3i are 
chosen positive and equal to /3 > 0, except for a fraction 
K of the bonds where the couplings are (— /3). In the 
zero temperature limit /? — > oo, the anti- ferromagnetic 
couplings induce frustration. 

It was first observed by Nishimori 0, 0| that bond 
distributions with vastly different values for k may still 
share the same thermodynamical properties. This is due 
to a Z2 gauge symmetry, which becomes transparent if 
we introduce link variables = sign{Pi). With this 
notation, the thermal energy is given by 



E{T) = 



•.= {xy) 



(2) 



The partition function, eq. and observables such as 
the thermal energy, eq. Q, are invariant under the fol- 
lowing change of bonds and spin variables: 

(j^'^{x) = Vl{x) a{x) 

uf^y) = n{x) U(xy) ^-\y) , (3) 

where the gauge transformation takes values in Z2, 
il{x) — ±1. In order to characterize the frustration of 
the model, we introduce the plaquette variable, P{p) = 
YieepUi, constructed from the given bond background. 
This definition is borrowed from lattice gauge theory, 
where a non-trivial value P{p) — —1 indicates that a 
Z2 vortex intersects the plaquette p. Notice that the 
variable P{p) is invariant under the Z2 gauge transfor- 
mation and the distribution of vortices is thus the 
proper measure to quantify the frustration of the model. 

With a given vortex content, there is still a large 
number of gauge equivalent bond distributions which 
share the same physical properties. In particle physics, 
the bond distribution with the minimal number of anti- 
ferromagnetic bonds is known as Landau gauge, 
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For this choice of bond distribution, the ground state is 
always uniform, 

a^^ — ffy^ = const. 'ix,y , (Landau gauge) . (5) 

To see this, we can use eqs. (O and to express the 

energy of a given spin configuration {cr} in a Landau 
gauge background as 

E[a] ^ - ^1'' - - E ^T""^ 

e={xy) e={xy) 

> - ^ ijf ^ = - ^ a^^ Uf"- a^^ = E[a^^] 



where we used the maximum condition eq. Q for the in- 
equality and the definition eq. (jSJ for the uniform ground 
state. 

Let us consider a particular distribution of anti- 
ferromagnetic bonds on the lattice. The ground state 
energy can be obtained as follows: (i) Calculate the po- 
sition of the vortices on the lattice; (ii) construct the 
minimal number Na of anti-ferromagnetic links which 
arc compatible with the given vortices, (iii) obtain the 
gauge transformation which casts the original bond 
distribution to the one obtained in step (ii). As shown 
above, the gound state in the bond distribution (ii) is 
uniform and its energy can be read off directly. 



Eo ^ Na - 2Ni, 



(6) 



where — dL"^ is the total number of bonds on the lat- 
tice. Transforming back to the initial bond distribution, 
the energy Eq remains unchanged while the ground state 
becomes (t° = The practical difficulty in this algo- 

rithm lies in step (ii) which is NP hard except for the case 
d = 2, where Edmond's algorithm 16, 17j provides an ef- 
ficient solution in polynomial time jlj. Figure Q] shows 
a random distribution of = 150 vortices and the cor- 
responding minimal number of frustrated bonds, which 
were obtained with Edmond's algorithm. Also shown is 
the ground state energy as function of the vortex den- 
sity p = Uy/L'^. The data comprise an average over 100 
vortex distributions for each value of p. 

In the following, we will construct a cluster algorithm 
which strongly reduces auto-correlations near criticality 
and which, in addition, is gauge invariant. For this pur- 
pose, we follow the derivation of the Swendsen-Wang 
cluster algorithm but firstly divide the bonds into 
those with ferromagnetic and anti-ferromagnetic cou- 
plings, respectively. For the ferromagnetic bonds, we 
make use of the identity 

exp{l3ax<Ty} = e^ [{I - q) + qS^^^y] , (7) 

with q = 1 — e^"^^ . A similar expression holds for the 
anti-ferromagnetic bonds, 

exp{-/3a,CTj,} = e'5 [(1 - g) + g (1 - . (8) 

The partition function can thus be written as 

Z = e^'^Y. n{^(^^) +9 '^-..J (9) 

+ (?(-/?,) {{\-q) + g(l-J,^,j]} , 

where is the total number of bonds on the lattice. 
The bond activation variable n £ {0, 1} is introduced via 
the identity 
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FIG. 1: Left panel: — 150 randomly distributed vortices (red squares) on a 30 x 30 lattice and the corresponding minimal 
matching with anti- ferromagnetic bonds (black lines) . Ground state energy per link: Eo — —0.837778(1). Right panel: ground 
state energy as function of the density of frustration p. 
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FIG. 2: Left panel: Thermal energy as a function of /3 for the configuration depicted in figure Right: Specific heat for a 
30 X 30 lattice with Uv randomly chosen vortices. 



The partition function is now recast to 



where 



(11) 



i={xy) 



P{a, n) =^ U \ 9{Ui) \Sn, o (1 - g) + <5„, i g 6, 



For the cluster update algorithm we perform subsequent 
updates of the bond variables ni and the spin variables 
ax ■ Let us first discuss the update and consider a spe- 
cific bond i = {xy). If € is ferromagnetic, Ug, = +\, the 
corresponding bond variable ni is set to zero if the spins 
attached to (. are anti-parallel, ax 7^ ctj/; if ax = ay we 
set ng = 1 with probability q. For the anti-ferromagnetic 
bonds we reverse this procedure: If the spins attached 
to £ are parallel, we set ni = 0, and otherwise we set 
ng = 1 with probability q. All spins connected by ac- 
tivated bonds are said to be part of one cluster. We 
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call this algorithm a "hybrid cluster algorithm" since the 
spins in one cluster generically have mixed orientations. 
Notice that only parallel spins connected by a ferromag- 
netic bond and anti-parallel spins connected by an anti- 
ferromagnetic bond can be part of the same cluster. 

Next we turn to the spin update: In a Swendsen-Wang 
type ^3 update, we would randomly choose i = ±1 and 
assign z to a randomly chosen spin ax- By definition, all 
spins in a^^s cluster can be connected to ax by a con- 
tinuous path of activated bonds. There may in fact be 
several such paths which involve a different number of 
anti-ferromagnetic bonds. However, if a particular con- 
nection of ax and ay on the cluster involves an even (odd) 
number of antiferromagnetic links, then every other clus- 
ter connection of ax and ay will also involve an even 
(odd) number. If our target spin ax is assigned the value 
I, then all cluster spins ay with even connections to x 
must be assigned the same value i, while all ay with odd 
connections must be assigned (— i), in order to avoid a 
configuration with zero probabilistic weight. In practice, 
we used a Wolff type variant Rather than growing 
all clusters on the lattice and flipping the spins with 50% 
probability, we pick a target spin, grow the correspond- 
ing cluster according to the rules above and the flip the 
entire cluster 

The only non-trivial part in the proof of the above 
algorithm is to show ergodicity, i.e. the fact that any 
spin configuration can be generated with a non- vanishing 
probability. In order to verify this, it is sufficient to show 
that the change of a single spin occurs with non- vanishing 
probability. This is possible e.g. if the algorithm identifies 
this single spin as a one-spin cluster, for which there is 
always a non-zero probability. 

Let us finally demonstrate that the algorithm is indeed 
gauge invariant. To this end, we assume that a partic- 
ular spin axg was selected by the algorithm to be part 
of a cluster. A local gauge transformation, ri(xo) = — 1, 



i}{x ^ Xo) = 1, changes axg — > —<^xo- At the same time, 
however, all bonds which are attached to xq also change 
their sign. From the cluster rules above, it can be seen 
that the spin ax„ would be part of the same cluster even 
after a local gauge transformation. Every gauge transfor- 
mation can be composed as a sequence of local one-spin 
transformations of the above type. Hence, the cluster 
growing prescription does not depend on the gauge. 

In order to demonstrate that the algorithm yields the 
same results after a gauge transformation, we have cal- 
culated the thermal energy as function of P for the spin- 
glass shown in figure ^ In Landau gauge, there are 146 
frustrated bonds to represent the riy = 150 vortices im- 
plying kl = 146/1800 « 0.081. After a random gauge 
transformation, this number increased to 888 frustrated 
links which roughly corresponds to a frustration density 
K « 0.5. Indeed, the algorithm provided the same result 
for these physically equivalent systems with vastly differ- 
ent K, cf. figure|2 We have also verified that hot and cold 
starts yield the same results in both cases. In addition, 
the dashed line in figure |21 indicates the ground state en- 
ergy of the spin glass, which is the energy in the limit 
f3 oo. As a first application, we calculated the specific 
heat for a particular spin-glass with randomly dis- 
tributed vortices. The result is shown in the right panel 
of figure El As expected, the (pseudo-)critical point is 
shifted to larger values of /3 if the number of defects 
is increased. At the same time, the peak broadens. 

In conclusions, we have stressed the importance of a 
gauge invariant classification of frustration. For the 2d 
frustrated Ising model, the exact ground state energy was 
discussed as a function of the density of gauge invariant 
vortices. Finally, a gauge invariant cluster update algo- 
rithm was developed which allows for efficient computer 
simulations near criticality. 
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